Quasiparticle and Optical Properties of Rutile and Anatase Ti02 
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Quasiparticle excitation energies and optical properties of Ti02 in the rutile and anatase struc- 
tures are calculated using many-body perturbation theory methods. Calculations are performed for 
a frozen crystal lattice; electron-phonon coupling is not explicitly considered. In the GW method, 
several approximations are compared and it is found that inclusion of the full frequency dependence 
as well as explicit treatment of the Ti semicore states are essential for accurate calculation of the 
quasiparticle energy band gap. The calculated quasiparticle energies are in good agreement with 
available photoemission and inverse photoemission experiments. The results of the GW calculations, 
together with the calculated static screened Coulomb interaction, are utilized in the Bethe-Salpeter 
equation to calculate the dielectric function e2{uj) for both the rutile and anatase structures. The 
results are in good agreement with experimental observations, particularly the onset of the main 
absorption features around 4 eV. For comparison to low temperature optical absorption measure- 
ments that resolve individual excitonic transitions in rutile, the low-lying discrete excitonic energy 
levels are calculated with electronic screening only. The lowest energy exciton found in the energy 
gap of rutile has a binding energy of 0.13 eV. In agreement with experiment, it is not dipole allowed, 
but the calculated exciton energy exceeds that measured in absorption experiments by about 0.22 
eV and the scale of the exciton binding energy is also too large. The quasiparticle energy alignment 
of rutile is calculated for non-polar (110) surfaces. In the GW approximation, the valence band 
maximum is 7.8 eV below the vacuum level, showing a small shift from density functional theory 
results. 

PACS numbers: 71.20.Nr 71.15.Qe 71.35.Cc 
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I. INTRODUCTION 

Even after half a century of researcbii^, investigation 
of the fundamental properties of titanium dioxide (Ti02) 
crystal phases remains important and fruitful, in part due 
to the role they have in concepts to effectively utilize solar 
energy. For example, Ti02 structures form the photoac- 
tive component in heterogeneous photo-catalysts which, 
by absorbing energy from the sunlight, degrade environ- 
mentally hazard materials^''* and split water into H2 and 
O^^. Scintered anatase Ti02 nanoparticles provide the 
backbone for electron transport and the substrate for 
organic chromophores in the Gratzel photovoltaic solar 
cells^. In addition to that, Ti02 has been widely used in 
various areas from optical coatings to pigments^. Funda- 
mental to all of these applications are the relative align- 
ments of essential energy levels near the valence and con- 
duction band edges of Ti02 crystal phases and the cor- 
responding optical transition energies. If predictive com- 
putational methods are going to have impact on the un- 
derstanding and design of heterogeneous photocatalytic 
systems based on Ti02, we must first establish that these 
methods can predict the basic properties of the crystal 
phases, providing a coherent framework for all the exper- 
imental facts. 

Rutile and anatase are two common crystal structures 
in which Ti02 is found. In both phases, each Ti atom in 
the crystal is surrounded by a slightly distorted octahe- 
dron formed by six oxygen atoms. The distinct phases 
exhibit a different connection between the distorted octa- 
hedra (TiOe)- In the rutile phase each octahedron shares 
two edges with its neighbors, while in the anatase phase 



each octahedron shares four^. In the rutile form, the 
crystal has a simple-tetragonal structure^ with a = b = 
0.45936 nm and c = 0.29587 nm. The symmetry of the 
lattice is described by the space group P4:2/mnm with 
the only internal parameter u = 0.30479. In the anatase 
form the crystal structure is body-centered tetragonal^iiS 
and belongs to space group lAi/amd. The three sides 
of the conventional cell are a = b = 0.3784 nm and c 
= 0.9515 nm respectively. The internal parameter u is 
0.2082ii£. The measurements quoted were done at room 
temperature; the change in lattice parameters upon re- 
ducing the temperature to 15 K is less than 0.001 nm^i.. 

Most of the early first-principles calculations of the 
properties of Ti02 were based on the local density 
approximation (LDA) in a Density Functional Theory 
(DFT) based approach^^'^^. The crystal structures and 
ground-state properties were accurately reproducedi^^— . 
However, as has been more generally observed for semi- 
conductors and insulators, the energy gaps pertaining to 
optical properties were found to be too small. The min- 
imum energy gap in the LDA band structure underesti- 
mated the band gap observed in optical experiment a^^i^^ 
by about 40%i2,. Calculations based on Hartree-Fock 
theory have been performed, giving accurate structural 
properties for rutile and anatase, but with a minimum 
energy gap that exceeded 10 eV^°. A hybrid approach, 
admixing a fraction of the bare exchange from Hartree- 
Fock, also showed accurate structural properties for 
rutile^*', with a band gap that is closer to experiment 
(3.4 eV)^""^. A more economical approach, approximately 
accounting for explicit Coulomb interactions through a 
U parameter acting on the Ti 3d electrons in a DFT+U 
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approach, overestimated the lattice parameters while still 
showing a band gap that was smaller than experiment^. 

A direct approach to calculate electronic excitation en- 
ergies based on the Green's function approach of many- 
body perturbation theory (MBPT), specifically utilizing 
the GW approximation for the electron self energy^^, 
has proven to be relatively accurate for a broad array 
of materials^ir— . Several calculations have been re- 
ported for Ti02 with different implementations of the 
GW method^^"— . However, the band gap was signifi- 
cantly overestimated in all these reports. 

The frequency dependent macroscopic dielectric func- 
tion probed in optical measurements has been exten- 
sively studied within the framework of DFT for Ti02 
crystal phase o^^i^^'^^i'^^ . Using the independent particle 
approach of Ehrenreich and Cohen^'^ , the underestimate 
of the fundamental band gap immediately gives an er- 
ror in the optical threshold and the overall shape of the 
dielectric function calculated in this way was quite dif- 
ferent from the experimental measurements^ii^^. In the 
framework of MBPT, the GW-based results for the elec- 
tronic excitation energies are input to a direct treatment 
of neutral excitations through solution of the two-particle 
Bethe-Salpeter equation (BSE), an approach that has 
provided a satisfactory description of the optical prop- 
erties of a number of systems^li^. An application of the 
BSE approach for rutile and anatase Ti02 has only re- 
cently appeared in the literature^^. The shape of the 
spectra are in much better agreement with experiment 
as compared to the independent particle approach. 

In this article, we critically assess the application of 
MBPT methods to calculate the electronic and optical 
excitations for Ti02 in the rutile and anatase crystal 
phases. To treat the electronic excitation energies, we 
use the GW method without self consistency or inclu- 
sion of vertex corrections. Empirically, this approach is 
often relatively accurate, although a full understanding 
of cancellations between self consistency and vertex cor- 
rections remains an open area of researcb^^^— . Several 
different approximations to handle the frequency depen- 
dence of the screened Coulomb interaction in the GW 
method are compared. We find that use of plasmon pole 
models^ii^ results in significant overestimation of the 
band gap. When the full frequency dependence of the 
screening is included, together with explicit treatment of 
the Ti semicore states, the calculated electronic excita- 
tion spectrum is found to be in good agreement with pho- 
toemission and inverse photoemission spectra^i^. In- 
terestingly, the calculated fundamental gap (3.34 eV and 
3.56 eV for rutile and anatase respectively), is still larger 
than the measured minimum gap from optical absorp- 
tion (3.03 eV for rutile^^ and estimated to be 3.3 eV for 
anatasei^) . As a first step towards application to hetero- 
geneous photocatalytic systems, we discuss the alignment 
of the valence and conduction band edges at non-polar 
vacuum-solid (110) interface of Ti02. We find that the 
GW method implies only very small corrections relative 
to the LDA for the valence band position, with the most 



of the band gap error going to shift the conduction band 
edge upwards. 

To explore the role of electron-hole interactions and ex- 
citonic binding energy, we have used the results from the 
GW based calculations as input to the BSE approach. 
Similar to the recent results of Lawler et al.^'^ , the calcu- 
lated frequency dependent dielectric function accurately 
reproduces the main onset of absorption near 4 eV and 
gives a good account of the frequency dependence for 
both rutile and anatase. We also solve the BSE for the 
low-lying, bound exciton states for rutile. The deepest 
exciton binding energy is calculated to be about 0.13 
eV. The dipole-forbidden character of the lowest exciton 
agrees with low temperature measurementsi^, although 
the predicted exciton energy (3.25 eV) is still 0.22 eV 
larger than experiments. Also, the magnitude of the ex- 
citon binding energy is larger. While the discrepancy 
for the exciton energy could very reasonably be regarded 
as within the expected errors of the MBPT methods 
used here, it may suggest an important role for electron- 
phonon coupling in screening and in further renormaliz- 
ing the energy gap in Ti02. 

The rest of the article is organized as follows: In Sec. 
II, the methodologies used in the DFT, electronic excita- 
tion and optical excitation calculations are briefly sum- 
marized, the key approximations are discussed and the 
numerical details are provided. In Sec. Ill, we present 
the main results for the electronic and optical excitations 
in rutile and anatase Ti02 and discuss them in compar- 
ison to available experiments. Finally, we conclude the 
article in Sec. IV with a short discussion, including the 
role of coupling to phonons. 



II. METHODOLOGY AND NUMERICAL 
DETAILS 

A. DFT Calculations 

The LDA eigenvalues and eigenvectors of Ti02 are 
calculated with a plane-wave basis set using norm- 
conserving pseudopotentials. Unless indicated otherwise, 
the LDA calculations are carried out using the ABINIT 
package^Silfi . In Ti02 the Ti is nominally ionized to 
[Ti^+l and the low lying conduction band states are of 
predominantly 3d character. As we show below, artifi- 
cially dividing the n=3 shell of Ti into frozen core (3s and 
3p) and valence {3d) contributions introduces a signifi- 
cant error to the energy band gap. The pseudopotential 
of Ti which includes semicore electrons is generated using 
the OPIUM package^ in the Troullier-Martins scheme^ 
with an initial configuration of (iVe)3s^3p^3(i°4s'^4p''. 
The outermost five orbitals are included and the cut- 
off radn (in Bohr) are 0.9, 0.9, 1.0, 0.9, and 0.9 re- 
spectively. Other pseudopotentials are taken from the 
ABINIT pseudopotential database generated using the 
FHI99PP package^^ 

In all calculations, the Perdew-Wang representation!^ 
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of Ceperly-Alder exchange-correlation potential^A is 
used. When including the Ti semicore states, a ki- 
netic energy cutoff of 200 Ry is used to ensure the con- 
vergence of the LDA results, as suggested by previous 
calculationsi^i^. To examine the accuracy of the pseu- 
dopotentials, we calculate the optimized lattice constants 
for rutile, finding a = 4.5484 A(4.5936 A), c/a =0.6414 
(0.64409) and u = 0.3040 (0.30479), agreeing with the 
experimental values noted in parentheses^ to the accu- 
racy generally expected for LDA calculations. We also 
compare our LDA calculations with results obtained us- 
ing the VASP package^^-^^ with the recommended pro- 
jector augmented-wave (PAW) pseudopotentials^s. The 
difference between the two LDA calculations is within 
0.5% for lattice parameters and less than 0.05 eV for 
bandgaps. In the GW and BSE calculations described 
below, the geometrical parameters of the unit cell for 
both rutile and anatase phases are taken from experi- 
mental measurements^"—!^. 



B. GW Method 

In MBPT, the evolution of the electrons in a material 
is described by the one-particle Green's function, with 
the effect of electron-electron interactions represented by 
an electron self energy operator. Well defined electronic 
excitations appear as peaks in the corresponding spec- 
tral function. Excitations with single particle character, 
namely quasiparticles, can be obtained as solutions of a 
Schrodinger-like equation^ 

{T + Ve,t+VH)^n^k {r) + 

J dr'j:{r,r';En,k)^n,k{r') = En,ki^n,k (1) 

where T is the kinetic energy, Vext is the external poten- 
tial, and Vh is the average Hartree potential. S is the 
self energy of the electrons and the indices refer to Bloch 
states n, k. It includes all the exchange-correlation effects 
contributed by surrounding electrons. Since E is gener- 
ally non-Hermitian, En^k is complex with the real part 
giving the quasiparticle energy and the imaginary part 
corresponding to the width of the quasiparticle peak in 
the spectral function, i.e. the quasiparticle lifetime. 

A practical approximation to calculate E has proven to 
be the so-called GW approximation of Hedin^S, in which 
the self energy I](r,r'; E) is formally written as 

E(r,r';i?) = 

— f dE'e-'^''^'G{r,r';E - E')W{r,r';E'), (2) 
27r J 



Here G is the Green's function of the electrons and W 
is the dynamically screened Coulomb interaction deter- 
mined by the inverse dielectric matrix e~^(r, r"; E), and 
(5+ is a positive infinitesimal time. The G and W in 
Eq. ^ refer to the fully interacting Green's function. 
However, in practice, using an initial LDA calculation to 
determine the screening through linear response calcula- 
tions (not including the exchange-correlation kernel) and 
to provide an initial, independent-particle Green's func- 
tion has often proven to be sufficiently accurate. There 
are specific examples where the LDA orbital character 
can be wrong, e.g. in some late 3d transition metal 
compounds^"—. However, in Ti02, the Ti 3d is almost 
empty and the valence band edge region is predominantly 
O 2p character with minimal admixture of Ti Sd. The 
Ti02 case should be similar to the vast majority of semi- 
conductors and insulators in this regard. Also, the LDA 
wavefunctions are sufficiently accurate that a first order 
estimate of the self energy correction to the LDA eigen- 
values is adequate. The quasiparticle energy correction 
AE„^k to a LDA orbital 0„.fc is obtained through a re- 
duced form of Eq. ([T]) as 

AEnM = ^n,fc (KM^iEnH^) - V^J'^\cl>n,k) , (3) 

where V^J^^ is the exchange-correlation potential and 
Zn_k is the renormalization factor of the orbital defined 
as Zn.k = (1 - dJ:/dE)-^\E=E^nA. 

The frequency dependence of the screened Coulomb in- 
teraction (W) can often be addressed using a generalized 
plasmon-pole (GPP) model^iii^, with substantial advan- 
tages in computational efficiency. The GPP models have 
proven to be relatively accurate for many semiconductors 
and insulators, including ionic crystals such as LiCl®^ and 
]y[gQ65 However, as discussed below, we find that use of 
the GPP leads to a gap that is substantially too large for 
Ti02. Several approaches to include the full frequency- 
dependent dielectric matrix have been implemented and 
described in the literature: (1) an analytical continua- 
tion method^iS, (2) a direct method which carries out 
the integration in Eq. ([2]) along the real axis^^"^°, and 
(3) a contour deformation (CD) method which deforms 
the integration in Eq. ^ along the imaginary axis^^. We 
adopt the contour deformation method to carry out the 
calculations, which is particularly efficient for evaluating 
self energy for states near the gap region. 

In the CD method, the correlation contribution 
T,'^ {r , r' ; E) of the self energy is written as the sum of 
two terms^iii^ 
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Wp{ry- \E - En^k] - ir]) MEf - En,k)e{En,k - E) - Q{Er,,k - Ef)Q{E - £;„,fe)]} , 



(4) 



where Wp = W — Vcoui, Ef is the Fermi energy, is a 
small damping amplitude and is a Heaviside function. 
The first term in Eq. comes from the integration 

along the imaginary axis. As W is now smooth along 
the imaginary axis, a sparse sampling of E is sufficient to 
converge the integration. The second term is the residual 
contribution of poles near the real axis. It is non-zero 
only while E > En^k > Ej ot E < En^k < Ef. For any 
E close to the Fermi surface, only Wp for \E — En.k\ ^ 
have non- vanished contributions to in the second 
term. This makes the calculation more computational 
efficient; is a smooth function of \E — En.k] around 0, 
due to the band gap, and only relatively low frequencies 
need to be sampled. 

We implement the contour deformation approach 
based on a private branch of the YAMBO package?-. 
The integration along the imaginary axis in Eq. Q is 
performed with a non-uniform mesh of N points accord- 
ing to 



El' = EotanC-^TT 



1,2 iV, 



(5) 



which maps the integration along the imaginary axis to 
an integration on the [0, 1) interval. The energy Eq pro- 
vides a scale for the overall density of the samples on 
the imaginary axis. Half of the mesh spans the energy 
scale from zero up to Eq while the other half sample the 
higher energies. For Ti02, a mesh of 50 points and an 
energy scale of 1 Ry were enough to keep the numerical 
error of the integration within 1 meV. Wp on the real axis 
is uniformly sampled with an energy increment of 0.1 eV 
and values between mesh points are linearly interpolated. 
The special case in Eq. ^ ior E ^ En^k must be han- 
dled with care. A consistent treatment, that avoids the 
apparent divergence and properly handles all the terms 
in Eq. (j4|), is to add a small positive energy to E — En.k 
(say S = 2.0 x lO"*" Ry) when necessary. A very similar 
contour deformation approach has been implemented in 
ABINIT~, and we have carefully compared the results 
for test cases. The ABINIT package has also been used 
for the GPP model calculations. 

For all the GW calculations, the energy cutoff is 60 Ry 
for the evaluation of the bare Coulomb exchange contri- 
bution , and 20 Ry for the correlation part Sc. A 
total of 160 bands are used for the calculation of both di- 
electric matrices and self-energies. An unshifted 4x4x6 
Monkhost-Pack (MP) mesh^** is used to sample the Bril- 
louin zone (BZ) of rutile, while for anatase an unshifted 
4x4x4 MP mesh is used. 

A test of the convergence with respect to the num- 
ber of bands included is shown in Fig. [T] for the fi- 
nal full-frequency approach with Ti semicore electrons 
treated explicitly as valence electrons in the pseudopo- 




Figure 1: Quasiparticle direct energy gap at F (a) and energy 
of VBM (b) for rutile as a function of the total number of 
bands kept in the full frequency dependent GW calculation. 
Squares represent calculations with Ti semi-core electrons in- 
cluded explicitly as valence electrons. The results are fitted 
using two different functional forms described in the text and 
displayed as solid lines that are indistinguishable in the figure. 
However, in (b) , they have different N ^ oo limits, displayed 
as dash-dot (l/N) and long dash {exp{—N/c)) lines in the 
figure respectively. 



tential. In order to characterize the fully converged val- 
ues, the data were fit with two different empirical forms, 
E{N) = Eo- b/N, and E{N) ^ Eo-b- exp{-N/c). We 
first check the validity of the fitting forms for the case 
of bulk silicon. The exponential form closely represents 
the band edge and band gap energies as a function of the 
number of included conduction bands, yielding extrapo- 
lated results in excellent agreement with those obtained 
via methods suggested by Bruneval and Gonze^^. For 
rutile, the fitting curves displayed in Fig. [T]are indistin- 
guishable, but predict slightly different — J> oo results 
for the absolute shift of the valence band edge, as indi- 
cated by the horizontal dashed lines in Fig. [IJb). In par- 
ticular, the fit for the quasiparticle energy gap indicates 
a converged quasiparticle energy of 3.37 eV for N=160. 
For the absolute shift of the valence band edge, the con- 
vergence is somewhat slower, with extrapolated values 
of -0.12 eV and -0.31 eV. This suggests that the valence 
band edge in the final results is 0.2 to 0.4 eV lower than 
the N=160 value. 



C. Bethe-Salpeter Equation and Optical Properties 

A detailed description of the BSE method is given 
in the literature^L^. We use the implementation in 
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the public branch of the YAMBO packageZ^. With the 
Tamm-Dancoff approximationl^ and the use of a static 
screened Coulomb interaction, the BSE is simplified to 
a generalized eigenvalue equation H^cv'c'^^ — EgA'^'^. 
The effective Hamiltonian H^c,v'c' has an explicit form of 



(6) 



where the quasiparticle energies (taken from the 
GW calculations) enter on the diagonal, Vvc,v'c' — 
{vc\Vcoui\v'c'), W^cyc' = {vv'\W{E^Q)\cc') and the 
indices w, c refer to the occupied valence and empty con- 
duction band states. For brevity, the explicit reference 
to Bloch wavevector k for each state is suppressed. The 
effective Hamiltonian in Eq. ([5]) is explicitly written for 
the spin singlet excitations which are directly probed in 
optical measurements. For the spin triplet excitations, 
the effective Hamiltonian is modified by dropping the so- 
called exchange term 2V^c,t)'c'- In terms of Eg and Ag'^, 
the macroscopic dielectric function eM(<^), including lo- 
cal field effects^liS^, is expressed as 



g->0 ^ — ' 



ljj - Eg 



IT] 



(7) 



In YAMBO, we use the option to evaluate the response 
function recursivel y^^i^^ . In order to study specific, low 
energy exciton states, we also directly diagonalize the 
generalized eigenvalue equation. 

Calculations of optical properties via BSE are more 
expensive computationally. For both phases, the static 
dielectric matrices are calculated with 80 bands and a 
damping coefficient of 0.1 eV. The dielectric function e(aj) 
is calculated on an 8 x 8 x 12 MP mesh for rutile and on 
an 12 X 12 X 12 MP mesh for anatase. For both cases, 
electron- hole pairs within 15 eV are taken in to build up 
the BSE kernel, in which the energy cutoff is 10 Ry for V 
and 3.5 Ry for W . To calculate the excitonic binding en- 
ergy of rutile, we restrict the basis set for the effective 
Hamiltonian to one conduction band and one valence 
band. The exciton binding energy converges relatively 
slowly with BZ sampling. The final results are reported 
based on a 12 x 12 x 18 MP mesh to sample the BZ. 
The energy cutoff is larger in this calculation for a more 
accurate representation of the BSE kernel, 14 Ry for V 
and 6 Ry for W. 



III. RESULTS 

A. Electronic Excitation Energies in Ti02 

The calculated electronic excitation energies in tita- 
nium oxides are found to be sensitive to technical factors 
in the GW calculations. We illustrate that here for the 
case of Ti02 in the rutile phase (Table |l|. First, ex- 
plicit, self consistent treatment of the semicore electrons 



(3s and 3p) on the Ti in the calculations for the solid 
affects the calculated energy gaps. As discussed in the 
literature, although the 3s and 3p levels are well sepa- 
rated from the id states energetically, there is significant 
spatial overlap^^. The effect for Ti02 is already evident 
at the LDA level, where the gap is reduced by more than 
0.3 eV upon explicit inclusion of the semicore electrons 
relative to freezing the semicore electrons in the pseu- 
dopotential. 

Depending on the approximation used to treat the elec- 
tron self energy in the GW method, the net influence 
of the semicore electrons varies. Focusing on the in- 
fluence of the semicore electrons for the most accurate, 
full- frequency method (FF), the change is greater than 
0.3 eV. A more detailed view of the contributions of the 
LDA exchange-correlation potential Vxc, the exchange S^, 
and the correlation Sc part of the self energy are also 
displayed Table H] The valence band maximum (VBM) 
orbitals are largely composed of oxygen p-states. They 
are spatially separated from the Ti semicore electrons, 
so the changes are relatively small. Explicit treatment of 
the semicore electrons as valence electrons only reduces 
Vxc by 0.24 eV. The reduction in Vxc is partially com- 
pensated by the decrease of Yi^ (0.08 eV) and Ec (0.04 - 
0.15 eV). Consequently the overall change of AEvbm is 
less than 0.1 eV. On the other hand, the conduction band 
minimum (CBM) orbital is almost purely of Ti 3d char- 
acter with substantial overlap to the semicore electrons. 
The expectation value of (Vxc) changes by about 9 eV. 
Of course, there is a corresponding, large change in the 
pseudopotential between these two cases. The bare ex- 
change term in the electron self energy operator changes 
by a similar amount. Their combined contribution to 
AEcBM is only about -0.34 eV. The changes in the cor- 
relation part of the self energy Ec of the CBM orbital 
is sensitive to the GW method. For the full-frequency 
method, Ec is decreased by only 0.1 eV when semicore 
electrons are included in the PP, so we flnd that most of 
the net effect on the band gap comes from the difference 
between LDA and bare exchange interactions with the 
semicore electrons. 

In Table H] the results of using different methods to ad- 
dress the frequency dependence of W are shown. These 
affect the correlation part of the electron self energy 
Ec and the renormalization factor Z. The results ob- 
tained with the full-frequency dependent dielectric func- 
tion, evaluated using the contour deformation method 
are the reference results, designated FF(CD) in the ta- 
ble. For comparison, results from two different general- 
ized plasmon pole models are shown. In the Hybertsen- 
Louie approach^!, designated GPP(HL), sum rules are 
applied to each individual dielectric matrix element to 
develop a plasmon-pole model for its frequency depen- 
dence. As shown in the lower part of Table HI when the 
semicore electrons of Ti were explicitly included in the 
pseudopotential, the calculated energy gap is almost 1 
eV too large. In the approach of von der Linden and 
Horsch^, designated GPP(VDLH), each dielectric ma- 
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Method SC VBM CBM 

inPP V,c Z Z LDA GW 



GPP(HL) 


n 


-19.97 


-24.59 


3.94 


0.82 


GPP(VDLH) 


n 


-19.97 


-24.59 


4.37 


0.82 


FF(CD) 


n 


-19.97 


-24.59 


4.50 


0.71 


GPP(HL) 


y 


-20.21 


-24.67 


4.09 


0.83 


GPP(VDLH) 


y 


-20.21 


-24.67 


4.41 


0.83 


FF(CD) 


y 


-20.21 


-24.67 


4.56 


0.70 



-11.48 


-3.57 


-5.69 


0.84 


2.08 


4.48 


-11.48 


-3.57 


-6.35 


0.85 


2.08 


3.60 


-11.47 


-3.57 


-5.78 


0.73 


2.08 


3.73 


-20.37 


-12.13 


-5.64 


0.85 


1.75 


4.27 


-20.37 


-12.13 


-5.99 


0.85 


1.75 


3.70 


-20.37 


-12.12 


-5.88 


0.72 


1.75 


3.38 



Table L Analysis of the valence band maximum (VBM), the conduction band minimum (CBM) and the direct energy gap Eg 
at the r point of the Brillouin zone calculated for rutile Ti02 using the GW method with several different approximations. 
The methods refer to two different generalized plasmon-pole (GPP) models and the full frequency-dependent (FF) approach 
described in the text. The second column indicates whether Ti semicore states are explicitly included as valence electrons. For 
the VBM and CBM, the expectation value is shown for the exchange-correlation potential Vxc in the LDA, the bare exchange 
self energy (S^,), the correlation part of the self energy (Ec) and the renormalization factor (Z). The band gap is shown in the 
LDA and for the GW method for each case. Energies are given in eV and the renormalization factor Z is dimensionless. 



trix is first transformed to the basis of eigenpotentials 
and then sum rules are applied to develop a plasmon 
pole model for each eigenpotential. This model also over- 
estimates the energy gap, but by a smaller amount. The 
renormahzation factor Z , in addition to entering the eval- 
uation of the quasiparticle energies in Eq. ([3]), indicates 
the degree of dynamical admixture of collective excita- 
tions into the quasiparticle. Larger deviations from unity 
indicate more admixture. Compared with the FF results, 
the value of Z indicates that the degree of admixture is 
substantially underestimated by both GPP models. 

We have also tested the GPP(HL) approximation for 
anatase Ti02 as well as two other titanates, SrTiOs and 
BaTiOa. In all three cases, the GPP(HL) approximation 
leads to minimum band gaps that are too large by 0.7-0.8 
eV. A similar deviation for the renormahzation constant, 
Z, is also observed. These results suggest that the quan- 
titative issues with the plasmon pole model extend more 
generally to titanates. In previous calculations of the loss 
function^ and the finite wavevector dynamical scatter- 
ing factor— for rutile Ti02 , substantial structure is seen 
in the frequency dependence, well beyond what could be 
easily accounted for by a single pole model. These ef- 
fects trace to an interplay between strong local field ef- 
fects and the Ti semicore p- to empty d-shell excitation. 
However, further analysis of the frequency dependence of 
the screening in Si and LiCl shows that deviations from 
a pole model for a range of frequencies around the plas- 
mon energy is not sufficient to predict the performance of 
the GPP model as it is used to evaluate the GW expres- 
sion. The dynamical screening at larger frequencies only 
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Table IL LDA and quasiparticle energy levels of rutile and 
anatase near the Fermi surface at selected k-points of high 
symmetry. Corresponding quasiparticle energy gaps are also 
displayed. The energy reference is taken to be the valence 
band maximum for both LDA and GW results. All energies 
are in eV. 



enters in an integrated fashion, resulting in substantial 
cancellations internally. In the case of titanates, we find 
that the full frequency dependence is essential for quan- 
titative results. Similar conclusions were drawn for the 
case of metallic Cu^. 
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The quasiparticle energy levels of rutile and anatase 
from the highest valence band and the lowest conduction 
band at selected high symmetry k-points are displayed in 
Table [TTl While the energy gap of rutile is found to be a 
direct gap (at F) in the LDA, our FF GW calculations 
indicate it as indirect (F — ^ R). However, the energy 
difference between the direct and indirect gap is small. 
The energy of CBM at F is only 0.04 eV higher than the 
energy at R. 

The calculated quasiparticle energies can be directly 
compared to spectroscopic measurements for electron re- 
moval or addition to the solid. The calculated value of 
the quasiparticle energy gap, 3.34 eV, agrees well with 
electron spectroscopy measurements, photo-emission and 
inverse photo-emission measurements^^i^. In Fig. [2l 
the density of states (DOS) of rutile derived from FF 
GW calculations is plotted together with the exper- 
imental spectra measured using x-ray photoemission 
and bremsstrahlung isochromat spectroscopy*^, which 
yielded a band gap of 3.3 ± 0.5 eV. Overall, the shape 
of the calculated DOS matches the shape of the exper- 
imental spectra, especially around the band gap. The 
experimental spectra measured using ultraviolet photoe- 
mission and inverse photoemission spectroscopy^ show 
very similar results with the inferred minimum energy 
gap about 0.2 eV smaller. 

The LDA calculations show that the top of the valence 
band of anatase lies in the A direction, somewhere about 
0.88 times of the distance from F to X. The energy at 
this k-point is 0.05 eV higher than the energy of the VBM 
at X. Subtracting the energy difference as a perturba- 
tion from the quasiparticle energy gap between X and F, 
which is 3.62 eV from Table HIl the quasiparticle energy 
gap of anatase is found to be 3.56 eV. The photoemis- 
sion data available for anatase show an overall occupied 
band with the oxygen p-bands^. That is similar to ru- 
tile. To our knowledge, there is no inverse photoemission 
data available for anatase, so the calculations can not be 
compared to a direct measurement of the quasiparticle 
energy gap. 



B. Optical Excitation Energies in Ti02 

More precise measurements of the minimum energy 
gap rely on optical absorption. This introduces the extra 
complication that the observed threshold for absorption 
will be altered by interactions between the photoexcited 
electron and hole, the formation of bound exciton states. 
For rutile, the BSE calculation shows a series of singlet 
bound excitonic states at F. The lowest two show s-state 
symmetry in the electron-hole envelope. They have a 
binding energy of 0.13 eV and 0.06 eV respectively. They 
are not dipole-allowed. The third and forth are degener- 
ate, with electron-hole envelope showing px,y symmetry, 
and have a binding energy of 26 meV. They are weakly 
dipole allowed for the electric vector perpendicular to 
the c-axis. Together with the calculated direct quasipar- 



ticle gap energy from above, we obtain the lowest energy 
singlet exciton at F with energy 3.25 eV and the first 
dipole allowed singlet exciton with energy of 3.35 eV. 
High resolution, low temperature optical absorption mea- 
surements for rutile resolve several separate features A 
very weak, but sharp exciton feature at 3.031 eV is identi- 
fied as the Is exciton which is electric quadrapole allowed. 
A stronger, but still relatively weak, dipole allowed 2p ex- 
citon feature starts at 3.034 eV. Finally, phonon-assisted 
features are also identified that correspond to an indirect 
gap of 3.049 eV. 

There are several important points of comparison. 
First, the calculated lowest exciton energy at F is about 
0.22 eV higher than measured. Broadly, this error is com- 
parable to those encountered when using the GW ap- 
proximation for other semiconductors^"—. However, it 
is important to be clear that the calculation is performed 
for a frozen lattice with no account given for electron- 
phonon interactions. In general, the electron-phonon in- 
teraction will reduce the zero-temperature quasiparticle 
gap^'^. Second, the present GW calculation gives the con- 
duction band minimum at the R point instead of being 
at F, as suggested by the optical measurements. The 
energy differences are small; the calculated conduction 
band at R is 0.04 eV lower than at F. In the analysis of 
the absorption data, the indirect gap is found to be 0.01 
eV higher than the direct ga.p^, albeit including what 
ever role excitonic effects may have. The difference be- 
tween theory and experiment is too subtle to be resolved 
in these calculations, particularly without the influence 
of electron-phonon interactions. Third, the symmetry of 
the excitonic states from the calculation agrees with the 
interpretation of the experiments. However, the scale of 
the excitonic effects that we calculate using the static 
dielectric matrix, and only including the electronic po- 
larizability, is substantially larger than suggested by the 
experiment. Our calculated long-wavelength dielectric 
constant ('^ 8) is slightly higher than the measured eoo 
(~ 7)^, but similar to previous calculations^. The lat- 
tice contribution is quite large, with eg = 111—. This 
again points to the importance of the electron-phonon 
interaction. 

The measured optical absorption in single crystal 
anatase at low temperature does not resolve any signifi- 
cant structure^. The energy dependence near the onset 
of absorption is consistent with an Urbach tail. Analysis 
of the temperature dependence leads to an estimate for 
the band gap for extended states of 3.42 eV^^. This ex- 
ceeds the measured optical threshold in rutile by about 
0.4 eV. Since the measured absorption edge in anatase 
is featureless, another way to characterize the absorp- 
tion edge and make comparison to rutile is to consider 
the energy at which the low temperature absorption is 
50 cm~^ for electric vector perpendicular to the c-axis. 
In rutile, this occurs at 3.04 eV while it occurs at 3.30 
eV in anatase^^. This suggests a more modest 0.26 eV 
difference between the minimum energy gap of anatase 
and rutile. Our calculated quasiparticle minimum energy 
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Figure 2: Density of states (DOS) of rutile derived from FF 
GW calculations plotted with photoemission and inverse pho- 
toemission spectra*^. The solid curve is the calculated DOS 
which is convoluted with a Gaussian function of (j= 0.5 eV, 
while circles are photo- and inverse photo-emission spectra. 



gap in anatase is 3.56 eV, modestly larger than the value 
deduced from the absorption measurements. The calcu- 
lated difference in gaps between anatase and rutile is 0.22 
eV, similar to the measured difference. 

So far, the calculated quasiparticle energies for Ti02 
have been found to agree well with electron spectroscopy, 
but the minimum energy gaps, including excitonic effects, 
are larger compared to the measured absorption thresh- 
old. Furthermore, the strength of the calculated excitoir 
binding energy is larger than implied by the interpre- 
tation of the optical spectrum near threshold. To get 
additional perspective, we calculate the macroscopic di- 
electric function over a broad energy range, includiirg the 
correlations induced by electroir-hole interactions. 

In Fig. [3l we show the imaginary part of the dielec- 
tric function of rutile for polarizations both perpendicu- 
lar and parallel to the tetragonal axis c. The solid curves 
are calculated from the BSE, while the dashed curves are 
derived from optical reflectivity measurements at room 
temperature^. For both polarizations, the theoretical 
spectra are close to experiment up to about 6 eV. In 
particular, the first strong peak at ^ 4 eV for both po- 
larizations is reproduced very well by the BSE results. 
Above 6 eV, the overall magnitude and prominence of 
the peaks in theoretical spectra are distinct from experi- 
ment. The 62(0;) of anatase is displayed in Fig. |31 where 
the experimental data'^'"' were measured at 100 K. Similar 
to the rutile case, the theoretical calculations capture the 
features around the onset of major absorption at 4 eV. 
For perpendicular polarization, the calculated oscillator 
strength is systematically too large starting at about 5 
eV. The calculated results are very similar to the previous 
calculation for rutile Ti02^. For anatase, the calculated 
peak heights near 4 eV appear less intense in their spec- 
tra, but this is largely due to their choice of a larger 
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Figure 3: e2(uj) of rutile from to 12 eV. Solid curves are 
theoretical calculations with BSE, and dashed curves are ex- 
perimental results^ obtained at room temperature. In (a) 
the direction of polarization is perpendicular to the tetrago- 
nal axis c, and in (b) the direction of polarization is parallel 
to axis c. 



damping parameter, as is evident from the broadening 
on the low energy side of their spectra. In particular, we 
have analyzed the integrated oscillator strength (i.e. the 
contribution to the f-sum rule) from the first peak in the 
anatase spectrum for parallel polarization. We find that 
our oscillator strength is essentially the same as theirs, 
but that both calculations show more oscillator streirgth 
than is found in the experimental spectra by about 30%. 



The systematic overestimation of oscillator strength 
at higher energy appears to be a more general issue. 
For example, a recent BSE study for several alkaline 
earth metal monoxides shows some similar excess oscilla- 
tor strength at higher photon energies^^. This may well 
trace to more fundamental assumptions in the method- 
ology. Two key issues are the use of the Tamm-Dancoff 
approximation, which has been identified recently as the 
main source of errors in the calculations of a confined 
system^^, and the assumption of a statically screened 
Coulomb interaction^iSi. Also, as noted by Lawler and 
coworkers, the f-sum rule for the oscillator streirgth con- 
verges very slowly in the titanates and the experimental 
analysis that relies on Kramers-Kronig transformations 
may have systematic errors as well~. 
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Figure 4: e2{uj) of anatase from to 12 eV. Solid curves are 
theoretical calculations from BSE, and dashed curves are ex- 
perimental results obtained at lOOK— . In (a) the direction of 
polarization is perpendicular to the tetragonal axis c, and in 
(b) the direction of polarization is parallel to axis c. 



External Energy Level Alignment: TiO2(110) 
Surface 



The energetic position of the quasiparticle VBM with 
respect to the vacuum level in the rutile phase for the 
(110) surface was calculated in several steps. First at 
the LDA level, analogous to the determination of work 
function for metals, the electrostatic potential change be- 
tween the bulk like region and the vacuum region was 
determined^^"— . Then in a second step, the bulk VBM 
position relative to the electrostatic potential is deter- 
mined. Previous calculations of work functions based on 
metal slabs based on the DFT Fermi energy were in fairly 
good agreements with experimental measurement a^^'^^ . 
In the third step, we will apply the perturbative correc- 
tion to the LDA exchange-correlation potential to deter- 
mine the alignment in the GW approximatio n^^i^^ . 

For semiconductors and insulators the Fermi level is 
separated from the VBM, and the position of VBM varies 
as the surface structure changes. To calculate the ener- 
getic alignment of rutile with respect to the vacuum from 
the (110) surface, a stoichiometric slab of six layers of Ti 
atoms is cut from the bulk structure with geometrical 
parameters described in Sec. I. A vacuum space of the 
same thickness of the slab is used to buffer the two sur- 
faces of the slab. The surface is cleaved as a (1x1), un- 
reconstructed surface as observed in experiments^. The 
k-point mesh is 2 x 4 x 1. After the slab is fully relaxed 
with the atomic positions of the central two layers fixed. 



the electrostatic potential inside the slab with respect to 
the vacuum is measured as -12.65 eV. The energy differ- 
ence between the VBM (from the LDA calculation) and 
the electrostatic potential inside the material is 5.05 eV, 
which is determined in a separate bulk calculation to pre- 
vent the quantum size effectsS^. Accordingly, the LDA 
value for the VBM relative to vacuum is 7.60 eV. Previ- 
ous LDA and gradient corrected calculations for give 7.16 
eV (LDA, three Ti layer slab)^, 7.2 eV (PBE, 11 Ti layer 
slab)92 and 7.6 eV (PW91, 11 Ti layer slab)S2,. In light 
of variations at the 0.2 eV level with number of layers in 
the slab model-^^, the overall agreement is satisfactory. 

The GW correction for the VBM calculated with FF 
model and 160 bands is 0.07 eV. However, as noted in 
Sect. lie, extrapolation to full convergence with respect 
to the total number of bands will drive this 0.2 to 0.4 
eV lower. We therefore suggest a GW correction of -0.2 
eV, with about 0.1 eV uncertainty. The final prediction 
from GW for the VBM alignment to vacuum at the clean 
rutile Ti02 (110) surface is 7.8 eV. 

In order to deduce the VBM alignment form experi- 
ment, two results must be combined: (1) the work func- 
tion which fixes the Fermi energy relative to the vac- 
uum; (2) the position of the VBM relative to the Fermi 
energy. Experimental measurements of the work func- 
tion of rutile from the (110) surface^^iSl vary from 4.7 
eV to 5.8 eV depending on the structure of the surface 
which is strongly influenced by treatment (annealing, ex- 
posure to oxygen, etc.). In addition, the position of the 
VBM relative to the Fermi energy also depends on sur- 
face treatmentiS^. Therefore, it is crucial to compare with 
data in which both values are measured on the same sam- 
ple. To our knowledge, this is relatively rare. Based on 
UPS measurements with {hi'=21.2 eV), a work function 
of 5.2 eV and a relative VBM position of 3 eV were mea- 
sured for TiO2(110)^'*, which implies the VBM position 
relative to vacuum falls at 8.2 eV. Similar measurements 
for TiO2(100) yield a work function of 4.9 eV and a rela- 
tive VBM position of 3.1 eV respectively^^, yielding the 
VBM position at 8.0 eV. The difference in workfunction 
between these measurements is consistent with a separate 
Auger Microprobe study of facet dependence^l. 

Based on this limited data set, the GW based pre- 
diction for the VBM alignment is off by about 0.4 eV. 
There are at least 0.1 eV uncertainties in both the the- 
ory and the experiment. Since both the work function 
and position of Fermi level are sensitive to the surface 
properties, the deviation between the theoretical calcu- 
lation and experimental measurements may well reflect 
the complexity of the Ti02 surfaces. For example, re- 
cent studies suggest that the commonly employed strat- 
egy of cleaning followed by annealing in oxygen may not 
result in the ideal surface envisioned (e.g. with no oxy- 
gen related defects)^. In particular, the physical origin 
of the widely observed defect states around 1 eV below 
the CBM remains a point of vigorous discussio n^^i^^ . 
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IV. CONCLUDING REMARKS 

We have presented a numerically well converged 
MBPT study of the electronic and optical excitation en- 
ergies in rutile and anatase crystals of Ti02. The calcula- 
tions are carried out in the approximation of a frozen lat- 
tice, without consideration of electron-phonon coupling. 
In most respects, the agreement with experiment is well 
within the expected accuracy of this approach. In par- 
ticular, the calculated quasiparticle gap agrees with elec- 
tron spectroscopy measurements (photoemission and in- 
verse photoemission) , the change in the minimum gap 
between rutile and anatase crystal structures is repro- 
duced, and the main features of the optical spectrum 
agree with ellipsometry measurements. The qualitative 
features of the zone center bound excitons calculated for 
rutile agree well with low temperature absorption mea- 
surements. However, the scale of the exciton binding 
energy is larger than that estimated from experiment by 
about a factor of 10 and the calculated exciton energy is 
about 0.2 eV larger than measured. 

The key theoretical assumptions in our application of 
MBPT include evaluation of the electron self energy in 
the GW approximation without iterating to self con- 
sistency or considering vertex corrections. Self consis- 
tency would certainly increase the calculated energy gap 
through the reduction in the screening22iii. Recent re- 
sults for a set of other semiconductors and insulators 
shows that approximate inclusion of vertex corrections 
in screening leads to a partially compensating reduction 
of the calculated energy gap^. However, a fully con- 
sistent approach to include vertex corrections remains 
subject of current research in the field. Because the Ti 
3d electrons are almost completely ionized, Ti02 should 
not be subject to the sorts of systematic errors found 
in non-selfconsistent calculations for some late transition 
3d metal compounds^"— The accurate results found for 
the key optical transition energies contributing to the 
absorption (Figs. 3 and 4) support this, and contrast 
to the case of CU2O where non-selfconsistent calcula- 
tions showed substantial discrepencies^. In the solution 
of the BSE, the calculated static (electronic) dielectric 
matrix has been used and the equations were simplified 
through the Tamm-Dancoff approximation. These ap- 
proximations are part of the standard MBPT treatment 
of optical spectra and the low energy excitons are ex- 
pected to be treated well21i2&. However, they may affect 
higher energy features in the spectrai^i. 

On physical grounds, we suggest that the most sig- 
nificant open issue concerns the role of electron-phonon 
coupling. In general, the electron-phonon self energy will 
both lead to a smaller zero-temperature quasiparticle gap 
and make a significant contribution to the temperature 
dependence of the energy gap*°. As noted in the text, 
there is a very large difference between the electronic 
dielectric constant and the low frequency dielectric 
constant including lattice polarization eg for Ti02. This 
suggests a relatively strong electron-phonon interaction 



and there is a long standing debate over the polaronic 
character of charge excitations in Ti02^. 

Based on the usual form of the Frohlich interaction, 
the dimensionless coupling constant characterizing elec- 
tronic coupling to the most important polar optic mode 
for rutile (with mode energy about 0.1 eV) is given by 
a = l.Qy^rrib/me where rri}, is the bare band mass and 
nie is the free electron mass. Using our DFT band dis- 
persions to have estimates, the electron (hole) band mass 
is about 0.6TOe (l.Sme) along x or y and l.Gmg (3.1me) 
along z. This suggests coupling constants of a ~ 1 — 2 
for electrons and a 2 — 3 for holes which would fall 
in the weak to intermediate coupling regime. Using the 
usual weak coupling expression, the electron and hole 
renormalization would be 0.1 — 0.2 eV and 0.2 — 0.3 eV 
respectively, both of which act to reduce the quasipar- 
ticle energy gap. For anatase a ^ ^.6^/mhJrn^, essen- 
tially the same as rutile, based on the mode energy and 
dielectric constants^^. The electron (hole) band mass 
is about 0.4me (l.Sme) along x or y and 3.9me(1.0TOe) 
along z, suggesting slightly different coupling constants 
of a '--^ 1 — 3 for electrons and a '--^ 2 for holes with cor- 
responding (weak coupling) electron and hole renormal- 
ization of 0.1-0.3 eV and 0.2 eV respectively. In weak 
coupling, the electron-phonon self energy is added to 
the results obtained here based on the GW approxima- 
tion. For the analysis of the optical absorption edge, 
a more detailed calculation is required because in the 
exciton-phonon coupling, the exciton is neutral and the 
electron and hole distortions of the lattice will partially 
canceliSLifl^. Taken together, if the large polaron regime 
is physically correct, these rough estimates suggest that 
the effect of the electron-phonon coupling could account 
for some of the differences between the present GW/BSE 
results for the frozen lattice and experiment. Firmer 
conclusions require a more extensive set of calculations, 
beyond the scope of this article22ii2^"i£^. In particu- 
lar, it may be that a more complete consideration of 
self-consistency and vertex corrections in the electron- 
electron contribution to the electron self energy will need 
to be combined with an analysis of the electron-phonon 
contribution. The two contributions should be treated 
in a fully consistent theory. In more empirical terms, 
an overstimate of the band gap based on selfconsistent 
treatment of the electron-electron interactions alone may 
be compensated by the electron-phonon contributions. 

It may well be the case that the large polaron regime 
is not applicable for Ti02. A recent THz spectroscopy 
study of rutile gave a direct measurement of the elec- 
tron scattering rate^°^. This data was analyzed with a 
Frohlich form for the electron-phonon interaction, but re- 
garding the coupling constant as a free parameter. Using 
the Feynman approachiSZiiS^ to handle intermediate to 
strong coupling, the analysis showed coupling constants 
for electrons a ~ 4 — 6 depending on field orientationiS^. 
The inferred electron mobilities were consistent with ear- 
lier electron transport measurement&i2^. These values 
suggest a substantially larger value for the electron self 
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energy of 0.4 — 0.7 eV-iS^. An older estimate based on 
small polaron theory also suggested 0.7 eV-i^^. A recent 
DFT+U based study suggested that an excess electron in 
rutile is indeed self trapped^^-*^. Although the binding en- 
ergy was not presented, the barrier for polaron hopping 
was estimated to be 0.3 eV. A similar study for an excess 
hole in rutile suggested barriers of 0.5 - 0.6 eV-i^S,. Taken 
together, if the small polaron regime is found to be phys- 
ically relevant, then the quasiparticle and excitonic ener- 
gies will need to be fully reanalyzed. For strong electron- 
phonon coupling a perturbative approach to combine the 
electron-electron and electron-phonon self energies is no 
longer justified. Furthermore, the electron-phonon cou- 
pling enters into the spectroscopic measurements in dis- 
tinct ways. The (inverse) photoemission and optical ab- 
sorption would each need to be properly analyzed. 
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